% hetero-nuclear PHIP program for 3 spins
format long e

% evaluating the energy levels of the 3 spin system
ai=sqrt(-1);

sx=[0 0.5; 0.5 0];                  % matrix of S_x for spin 1/2
sy=[0 -0.5*ai; 0.5*ai 0];           % matrix of S_y for spin 1/2
sz=[0.5 0; 0 -0.5];                 % matrix of S_z for spin 1/2
od=[1 0; 0 1];                      % 2*2 unity matrix

% setting {x,y,z} components of spin operators of all spins: 1024*1024
% matrices
I1z=zeros(1024,1024);
I2z=zeros(1024,1024);
I3z=zeros(1024,1024);
I4z=zeros(1024,1024);
I5z=zeros(1024,1024);
I6z=zeros(1024,1024);
I7z=zeros(1024,1024);
I8z=zeros(1024,1024);
I9z=zeros(1024,1024);
I10z=zeros(1024,1024);

I1x=zeros(1024,1024);
I2x=zeros(1024,1024);
I3x=zeros(1024,1024);
I4x=zeros(1024,1024);
I5x=zeros(1024,1024);
I6x=zeros(1024,1024);
I7x=zeros(1024,1024);
I8x=zeros(1024,1024);
I9x=zeros(1024,1024);
I10x=zeros(1024,1024);

I1y=zeros(1024,1024);
I2y=zeros(1024,1024);
I3y=zeros(1024,1024);
I4y=zeros(1024,1024);
I5y=zeros(1024,1024);
I6y=zeros(1024,1024);
I7y=zeros(1024,1024);
I8y=zeros(1024,1024);
I9y=zeros(1024,1024);
I10y=zeros(1024,1024);

% setting matrices of scalar products of all spins of interest
pro13=zeros(1024,1024);
pro14=zeros(1024,1024);
pro15=zeros(1024,1024);
pro23=zeros(1024,1024);
pro24=zeros(1024,1024);
pro25=zeros(1024,1024);

pro16=zeros(1024,1024);
pro17=zeros(1024,1024);
pro26=zeros(1024,1024);
pro27=zeros(1024,1024);
pro67=zeros(1024,1024);

pro68=zeros(1024,1024);
pro79=zeros(1024,1024);
pro69=zeros(1024,1024);
pro78=zeros(1024,1024);

pro610=zeros(1024,1024);
pro710=zeros(1024,1024);
pro810=zeros(1024,1024);
pro910=zeros(1024,1024);


Ix=zeros(1024,1024);
Iy=zeros(1024,1024);
Iz=zeros(1024,1024);

I1z=kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I2z=kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I3z=kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I4z=kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I5z=kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,od)))))))));
I6z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,od)))))))));
I7z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,od)))))))));
I8z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,od)))))))));
I9z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,od)))))))));
I10z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sz)))))))));

I1x=kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I2x=kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I3x=kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I4x=kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I5x=kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,od)))))))));
I6x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,od)))))))));
I7x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,od)))))))));
I8x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,od)))))))));
I9x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,od)))))))));
I10x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sx)))))))));

I1y=kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I2y=kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I3y=kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I4y=kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
I5y=kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,od)))))))));
I6y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,od)))))))));
I7y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,od)))))))));
I8y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,od)))))))));
I9y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,od)))))))));
I10y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sy)))))))));


p13=I1x*I3x+I1y*I3y+I1z*I3z;
p14=I1x*I4x+I1y*I4y+I1z*I4z;
p15=I1x*I5x+I1y*I5y+I1z*I5z;
p23=I2x*I3x+I2y*I3y+I2z*I3z;
p24=I2x*I4x+I2y*I4y+I2z*I4z;
p25=I2x*I5x+I2y*I5y+I2z*I5z;

p16=I1x*I6x+I1y*I6y+I1z*I6z;
p17=I1x*I7x+I1y*I7y+I1z*I7z;
p26=I2x*I6x+I2y*I6y+I2z*I6z;
p27=I2x*I7x+I2y*I7y+I2z*I7z;
p67=I6x*I7x+I6y*I7y+I6z*I7z;

p68=I6x*I8x+I6y*I8y+I6z*I8z;
p79=I7x*I9x+I7y*I9y+I7z*I9z;
p69=I6x*I9x+I6y*I9y+I6z*I9z;
p78=I7x*I8x+I7y*I8y+I7z*I8z;

p610=I6x*I10x+I6y*I10y+I6z*I10z;
p710=I7x*I10x+I7y*I10y+I7z*I10z;
p810=I8x*I10x+I8y*I10y+I8z*I10z;
p910=I9x*I10x+I9y*I10y+I9z*I10z;

p18=I1x*I8x+I1y*I8y+I1z*I8z;
p28=I2x*I8x+I2y*I8y+I2z*I8z;
p19=I1x*I9x+I1y*I9y+I1z*I9z;
p29=I2x*I9x+I2y*I9y+I2z*I9z;

p110=I1x*I10x+I1y*I10y+I1z*I10z;
p210=I2x*I10x+I2y*I10y+I2z*I10z;



Iz=I1z+I2z+I3z+I4z+I5z+I6z+I7z+I8z+I9z+I10z;
Ix=I1x+I2x+I3x+I4x+I5x+I6x+I7x+I8x+I9x+I10x;
Iy=I1y+I2y+I3y+I4y+I5y+I6y+I7y+I8y+I9y+I10y;



% J-couplings in Hz
J13=7.55;
J14=7.55;
J15=7.55;
J23=7.55;
J24=7.55;
J25=7.55;

J16=-0.6;
J17=-0.6;
J26=-0.6;
J27=-0.6;
J67=1.5;

J68=7.5;
J79=7.5;
J69=0.8;
J78=0.8;

J610=1.3;
J710=1.3;
J810=7.5;
J910=7.5;

J18=0.29;
J28=0.29;
J19=0.29;
J29=0.29;

J110=-0.45;
J210=-0.45;




tp=100.0;   % preparation time of PHIP: should be long that all coherences are washed out
te=0.0;     % free evolution time of PHIP at the preparation field, should be usually set to zero

% all chemical shifts in ppm in the sigma-scale
ch1=-2.332;
ch2=-2.332;
ch3=-1.053;
ch4=-1.053;
ch5=-1.053;

ch6=-7.207;
ch7=-7.207;
ch8=-7.265;
ch9=-7.265;
ch10=-7.155;

% Magnetic fields
% FV: time and number of points:
T2=2.0; % T_2 value, which gives the NMR linewidth 

[BotT] = textread('Bott01mT.DAT');  % B(t) profile of field switching
Nfv=2995;                           % number of points in the profile
        
B=10000*BotT(1,2);                  % starting (preparation) field
Bdet=70400;                         % final (detection) field

H = zeros(1024,1024);               % Hamiltonian: 1024/1024 matrix
    rho = zeros(1024,1024);         
    rho0 = zeros(1024,1024);
    rhot = zeros(1024,1024);
    
    rhos=zeros(4,4);
    rhos(2,2)=1/2;
    rhos(3,3)=1/2;
    rhos(2,3)=-1/2;
    rhos(3,2)=-1/2;
    

    H=B*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z+ch7*I7z+ch8*I8z+ch9*I9z+ch10*I10z)/7.04e4;
    H=H+J13*p13+J14*p14+J15*p15+J23*p23+J24*p24+J25*p25+J16*p16+J17*p17+J26*p26+J27*p27+J67*p67;
    H=H+J68*p68+J79*p79+J69*p69+J78*p78+J610*p610+J710*p710+J810*p810+J910*p910;
    H=H+J18*p18+J28*p28+J19*p19+J29*p29+J110*p110+J210*p210;
    
% setting starting density matrix: 1024*1024; spins 1 and 2 are ALWAYS
% coming from para-hydrogen
rho0=kron(od,kron(rhos,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od))))))));
        rho0=rho0/256;
        

% calculating the transition matrix 


    

% calculating the eigenvenctors
evec=zeros(1024,1024);
evec1=zeros(1024,1024);
energ=zeros(1024,1024);
HH=zeros(1024,1024);
[evec,energ]=eig(H);    % solving the eigen-problem for the Hamiltonian




evec;
evec1=inv(evec);  
        HH=evec1*H*evec;
        
        rho0=evec1*rho0*evec;   % transforming the DM to the eigen-basis

%check diagonalization!!!        
HH;

%DM evolution at polarization field: washing out the coherences during the
%preparation period

      for i=1:1024
          for j=1:1024
              
              ex=energ(i,i)-energ(j,j);
if ex==0
    rho(i,j)=rho0(i,j);
else
    rho(i,j)=rho0(i,j)*(exp(-2*pi*ai*ex*tp)-1)/(-2*pi*ai*ex*tp)*exp(-2*pi*ai*ex*te);
end              
          end
      end

rhot=rho;

rhot=evec*rhot*evec1;   % now going back to the Zeeman basis; new DM takes into account the evolution at the preparation field


% DM at polarization field is calculated
% now field variation
    Hi = zeros(1024,1024);


for ifv=1:Nfv-1

    Bint=10000*BotT(ifv,2);
    dt=BotT(ifv+1,1)-BotT(ifv,1);
    
    Hi=Bint*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z+ch7*I7z+ch8*I8z+ch9*I9z+ch10*I10z)/7.04e4;
    Hi=Hi+J13*p13+J14*p14+J15*p15+J23*p23+J24*p24+J25*p25+J16*p16+J17*p17+J26*p26+J27*p27+J67*p67;
    Hi=Hi+J68*p68+J79*p79+J69*p69+J78*p78+J610*p610+J710*p710+J810*p810+J910*p910;
    Hi=Hi+J18*p18+J28*p28+J19*p19+J29*p29+J110*p110+J210*p210;
    
    
   rhot=expm(-2*pi*ai*Hi*dt)*rhot*expm(2*pi*ai*Hi*dt);  % evolution during switching; Hamiltonian is Hi and time period is dt
    
    
    ifv
end

H0=zeros(1024,1024);    % Hamiltonian at the detection field; needed to calculate the NMR spectrum

    H0=Bdet*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z+ch7*I7z+ch8*I8z+ch9*I9z+ch10*I10z)/7.04e4;
    H0=H0+J13*p13+J14*p14+J15*p15+J23*p23+J24*p24+J25*p25+J16*p16+J17*p17+J26*p26+J27*p27+J67*p67;
    H0=H0+J68*p68+J79*p79+J69*p69+J78*p78+J610*p610+J710*p710+J810*p810+J910*p910;     
    H0=H0+J18*p18+J28*p28+J19*p19+J29*p29+J110*p110+J210*p210;
    
    
fi=pi/4;            % flip angle in NMR
    
    fx=zeros(1024,1024);
    fxx=zeros(1024,1024);
    
    fx=expm(ai*fi*Iy);
    fxx=expm(-ai*fi*Iy);
    
    [e0,v0]=eig(H0);
    e01=inv(e0);
    
    rhot=e01*rhot*e0;   % DM for PHIP-NMR
    rhoeb1=e01*Iz*e0;   % DM for equilibrium NMR
    fx=e01*fx*e0;
    fxx=e01*fxx*e0;
    magn=e01*(Ix+ai*Iy)*e0;
    
    AAA=zeros(1024,1024);
    BBB=zeros(1024,1024);
    
    % destroying all coeherences in the DM after the spin system arrives to
    % the detection field
                for i=1:1024
                for j=1:1024
               if (abs(i-j)>0)
                   rhot(i,j)=0;
                   rhoeb1(i,j)=0;
               end
                end
                end
    
ZA=fxx*rhot*fx;
ZB=fxx*rhoeb1*fx;
magnt=transpose(magn);

                
    for i=1:1024
        for j=1:1024
AAA(i,j)=AAA(i,j)+real(ZA(i,j)*magn(j,i));
BBB(i,j)=BBB(i,j)+real(ZB(i,j)*magn(j,i));
                end
    end

 
    
    Nomeg=10000;    % number of pints in the NMR spectrum
    omegF=-7.4;     % staring and final NMR frequency in the spectrum
    omegL=-0.8;
    
spe=zeros(Nomeg,3);
% now calculating NMR and PHIP-NMR
    
    for iw=1:Nomeg
    
            w=omegL-(iw-1)*(omegL-omegF)/(Nomeg-1);
            spe(iw,1)=w;
            w1=Bdet*w*3.0013e2/7.04e4;
for i=1:1024
for j=1:1024
    
    xx=abs(w1-v0(j,j)+v0(i,i));
    k2=1/T2;
    
   spe(iw,2)=spe(iw,2)+AAA(i,j)*k2/(xx*xx+k2*k2);
   
%   real(1.0/(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
      spe(iw,3)=spe(iw,3)+BBB(i,j)*k2/(xx*xx+k2*k2);
      
%      /(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
end
end

    end

    % plotting the result
plot(spe(:,1),real(spe(:,2)),spe(:,1),real(spe(:,3)));

% saving the result
save full_01mT.dat spe -ASCII